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Abstract 

The  inner  (singular)  integral  in  the  inverse  Radon  transform  for  paral¬ 
lel  beam  computerized  tomography  devices  can  be  integrated  analytical¬ 
ly  if  the  Radon  transform  considered  as  a  function  of  the  ray  position 
along  the  detector,  is  a  cubic  polynomial  spline.  Furthermore  by  using 
some  spline  identities,  large  terms  that  cancel  can  be  eliminated  ana¬ 
lytically  and  the  calculation  of  the  resulting  expression  for  the  inner 
integral  done  in  a  numerically  stable  fashion.  We  suggest  using  smoo¬ 
thing  splines  to  smooth  each  set  of  projection  data  and  by  so  doing  ob¬ 
tain  the  Radon  transform  in  the  above  spline  form.  The  resulting  analy¬ 
tic  expression  for  the  inner  integral  in  the  inverse  transform  is  then 
readily  evaluated,  and  the  outer  (periodic)  integral  is  replaced  by  a 
sum.  The  work  involved  to  obtain  the  inverse  transform  appears  to  be 
within  the  capability  of  existing  computing  equipment  for  typical  large 
data  sets.  In  this  regularized  transform  method  the  regularization  is 
controlled  by  the  smoothing  parameter  in  the  splines.  The  regularization 
is  directed  against  data  errors  and  not  to  prevent  unstable  numerical 
operations.  Strip  integral  as  well  as  line  i ntegral -data  can  be  hand¬ 
led  by  this  method.  The  method  is  shown  to  be  closely  related  to  the 
Tihonov  form  of  regularization. 

I ,  Introduction 

Consider  a  thin  "slice"  of  the  human  head.  In  modern  computerized 
tomography  (CT)  with  parallel  beam  geometry  the  equivalent  of  an  array 
of  2N  +  1  X-ray  beams  is  directed  through  the  slice  and  the  amount  of 
attenuation  of  each  beam  is  measured.  This  procedure  is  repeated  as 

the  array  is  rotated  through  M  positions,  s1 . sM,  about  the  head 

(see  Fig.l)  to  give  attenuation  factors  for  a  total  of  n*(2N+l)M  beams 
through  the  slice.  The  log  of  the  attenuation 
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factor  for  the  i th  beam  when  the  array  is  in  the  j th  position  is 
given  approximately  by 


Lijf  =  /  w..(x,y)f(x,y)dxdy 


(l.D 


where  f(x,y)  is  the  X-ray  density  of  the  head  slice  at  the  point 
(x,y)  and  w^  is  a  non-negative  weight  function  which  is  0  outside 
a  strip  surrounding  the  ij  th  beam  and  represents  the  non-uniform 
effective  distribution  of  the  beam  intensity  across  its  narrow  width. 
The  formula  makes  the  approximation  that  the  X-ray  attenuation  coeffi¬ 
cient  is  constant  over  the  spread  of  energies  present  in  the  (nearly) 
monochromatic  beam.  In  this  report  we  model  the  data  as 

2ij  =  Lijf+Eij  •  iB'N . N’ 

j  =  l  ,2 . M, 

where  the  c^.  are  independent  zero  mean  random  variables  with  approxi 
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mately  the  same  variance  which  model  counting  noise  and  any  other 
(hopefully  non-systematic)  errors  inherent  in  the  measuring  device 
and  the  approximations  being  made.  The  number  n(2N+l)M  of  data  points 
may  be  of  the  order  of  magnitude  of  10^. 

In  practice,  various  methods  are  used  to  estimate  f  from  the  data 
vector  z=(z_N  ^.....Zj^  M).  The  results  are  usually  presented  on  a 
video  display  with  different  values  of  the  estimate  of  f(x,y)  repre¬ 
sented  by  different  levels  on  a  gray  scale.  For  more  detailed  discus¬ 
sions  of  the  subject,  see,  for  example  Shepp  and  Logan  (1974),  Herman 
and  Naparstek  (1977). 

In  Section  2  we  review  briefly  the  estimation  of  a  function  f  by 
Tihonov  regularization  given  data  z  .*1  .f  ,  i*1 ,2  , . . . ,  n ,  where  the 
L.  are  arbitrary  continuous  linear  functionals  on  some  appropriate 
Hilbert  space.  This  approach  is  not  usual  in  human  head  and  body  CT 
because  of  the  apparent  numerical  difficulty  and  the  computational 
convenience  of  transform  methods.  See,  however  Natterer  (1980),  Artzy, 
Elfving  and  Herman  (1979).  For  our  purposes,  a  close  examination  of 
this  form  of  regularization  will  serve  to  clarify  the  resol uti on-noi - 
se  sensitivity  tradeoff  common  to  most  regularization  methods  for  dea¬ 
ling  with  discrete,  noisy  data.  The  method  is  highly  appealing  in 
many  mildly  ill  posed  problems  (as  is  the  CT  problem)  whenever  it  is 
feasible  to  implement  it. 

Most  modern  human  CT  devices  use  methods  for  estimating  f  based  on 
an  approximate  numerical  evaluation  of  a  regularized  inverse  Radon 
transform.  For  a  recent  description  of  one  such  algorithm,  see  Herman 
Naparstek  (1977),  Chang  and  Herman  (1978).  In  Section  3  we  present  a 
new  approach  for  the  approximate  numerical  integration  of  the  inverse 
Radon  transform  from  discrete,  noisy  data.  The  work  was  motivated  by 
a  study  of,  but  is  apparently  quite  different  from  the  method  described 
in  the  above  two  references.  It  is  in  fact  quite  close  to  the  Tihonov 
form  of  regularization  with  moment  discretization.  The  method  entails 
using  a  cubic  smoothing  spline  to  obtain  a  smooth  function  represen¬ 
ting  each  set  of  projection  data,  that  is,  each  set 
Zj  =  (z_N  j  • z  _  (  n_  i )  j  >  •  •  •  •  j)  where  j  is  fixed.  Then  the  inner 
(singular)  integral  in  the  Radon  transform  can  be  evaluated  analytically. 
After  using  certain  relations  between  the  coefficients  in  cubic  splines, 
one  obtains  a  computationally  stable  numerical  inversion  formula  which 
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appears  feasible  to  implement  with  data  sets  with  N  and  M  of  the 
order  order  of  10^*®. 

The  smoothing  parameter  in  the  cubic  smoothing  splines  controls  the 
resolution  noise  sensitivity  tradeoff.  The  suggested  approach  bypasses 
most  of  the  usual  discretization,  quadrature  and  aliasing  errors  common 
to  other  methods.  Unlike  smoothing  approaches  which  are,  at  least  in 
part,  directed  against  numerical  problems  connected  with  evaluating  a 
singular  integrand,  the  present  approach  directs  the  smoothing  against 
the  noisy  data  in  such  a  way  that  the  singularity  can  be  integrated 
out  analytical ly. 

In  Section  4  we  indicate  the  relationship  between  the  transform  method 
proposed,  and  Tihonov  regularization. 

2.  Tihonov  regularization 

Let  H  be  a  Hilbert  space  of  functions  defined  on  some  domain  n,  let 
fe H  and  suppose  one  observes 

zi  =  L.f+c.  ,  1-1,2 . n  (2.1) 

where  the  li  are  continuous  linear  functionals  on  H  ,  and  the  e^  are 
errors.  It  is  supposed  that  the  are  uncorrelated  zero  mean  random 

variables  with  common  variance. 

Having  chosen  w,  the  (Tihonov)  regularized  estimate  f  .  of  f  given 

n  i  a 

the  data  z-(z^ , . . .  ,z  )’  Is  the  solution  to  the  problem:  Find  feH  to 
minimize 

J-  .Zi  (Ltf-2  . )2  +  X||f||  2  .  (2.2) 

The  first  term  represents  the  "infidelity"  of  the  solution  to  the  data 

p 

and,  assuming  H  is  a  space  of  "smooth"  functions,  ||f  .||  represents 

n  j  a 

the  "roughness"  of  the  solution.  The  parameter  X  controls  this  tradeoff. 
Equivalently  X  controls  the  tradeoff  between  sensitivity  to  noise,  and  resolution. 

If  x  is  large  ||  f  x||  will  be  small,  and  the  solution  will  have  low 
resolution  but  the  sensitivity  to  noise  will  also  be  less,  since 
•»  n  o 

-  T  (L.f  ,-z<)  can  be  larger.  A  small  X  will  allow  || f  .  ||  to  be 

ni=!  1  n’*  1  n,A 

large  and  corrrespondingly  require  L.fn  x  to  match  the  data  better 
in  mean  square. 
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Since  the  are  continuous  linear  functionals  on  H  ,  there  exist 

representers  nj»...»nncN  such  that  L^f=<n^»f>  ,  where  <.  ,  .> 
is  the  inner  product  in  H.  Then  the  minimizer  of  (2.2)  can  be 
shown  to  be  given  by 

fntXm  <n1 . nn)(Q+nU)'1z  ,  (2.3) 

where  Q  is  the  n*n  gram  matrix  of  the  representers,  with  ijth 
entry  Qy 


Equivalently, 

f„,»  ■  K;<wnU>',!  •  <*•«) 

where  Kn:  K  -*En  is  defined  by  Knf*(L]f ... .  ,L  f) ,  and  K*  is  the 
adjoint  of  in  the  sense  that  (z ,Knf )=<K*z ' ,f>,  all  zeEn,feW 

where  (  ,  )  is  the  Euclidean  inner  product.  (KpK*  is  the  operator 
of  matrix  multiplication  by  Q)  Results  are  available  concerning  the 
convergence  of  fR  ^  when  x  is  chosen  appropriately  and  arc-  stated 
in  a  little  more  detail  in  Section  4. 

We  remark  that  if  H=t-2  then  K^KnKn+XI^"1  1  s  essentially  a  back 
projection  operator,  see  Natterer  (1980),  however  in  this  case  X  should 
be  thought  of  as  controlling  the  scale  or  dynamic  range  of  the  solution 
rather  than  its  smoothness,  and  it  is  then  not  very  important  parameter 
for  tumor  detection. 

We  make  some  remarks  on  choosing  X  and  the  space  H.  Natterer  (1980) 
has  suggested  that  for  computerized  tomography,  H  should  be  chosen  as 
the  space  Ha(S5), 

Ha(Jl)  *  ( f : // ( 1  + 1  E  |  ^  )a  |f(0|*d£  <  °°»  supp  fcjj) 

A 

where  f(^)  is  the  Fourier  transform  of  f  and  a  is  close  to  1/2. 
Ideally,  one  should  choose  H  so  that  it  "just"  contains  the  true  solu¬ 
tions.  If  one  looks  at  the  problem  in  "frequency  Space"  (see  Craven  and 
Wahha  (1979))  or  "eigenfunction  space"  (see  Wahba  (1979a)),  one  can  see 
that  the  regularized  estimate  fnX  can  be  thought  of  as  passing  the 
data  through  a  "low  pass  filter"  where  X  controls  the  half  power  point 
(or  "bandwidth")  of  the  filter  and  a  controls  the  "shape",  or  steepness 


A 
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of  the  "roll  off"  of  the  filter.  For  H  fixed  the  method  of  genera¬ 
lized  cross  validation  (GCV)  can  be  used  to  estimate  a  good  value  of 
X,  or  in  the  case  of  computerized  tomography,  to  obtain  good  starting 
values  for  human  “fine  tuning".  See  Wahba  (1979b)  and  references 
cited  there.  In  the  typical  tomography  problem  it  will  be  necessary 
to  utilize  the  special  structure  of  the  problem  and  possibly  to  do 
GCV  on  a  subset  of  the  data.  See  the  appendix. 

3.  A  novel  regularized  transform  method  using  smoothing  splines. 

Herman  and  Naparstek  (1977)  and  Chang  and  Herman  (1978)  have  recently 
studied  regularized  transform  methods  for  CT  reconstruction  for  a 
fan  beam  device.  In  this  section  we  suggest  a  new  numerical  approach 
to  the  regularized  inversion  of  the  Radon  transform  for  a  parallel  beam 
device.  A  similar  but  more  involved  analysis  can  be  carried  out  for 
the  fan  beam  inverse  transform  discussed  by  Herman  and  Naparstek  but 
we  do  not  do  it  here.  The  method  to  be  given  appears  to  have  the 
advantage  of  introducing  discretization  errors  and  quadrature  appro¬ 
ximations  relatively  late  in  the  numerical  procedure,  and,  intuitively, 
the  regularization  parameter  of  the  method  appears  to  affect  the  reso¬ 
lution  -  noise  sensitivity  tradeoff  in  an  appropriate  manner.  The  ncise 
suppresion  filtering  acts  directly  on  the  raw  data.  The  resulting 
smoothed  data  is  in  such  a  form  that  the  singular  integrand  is  evalua¬ 
ted  analytically,  and  large  terms  which  cancel  are  subtracted  analytically. 
Unstable  numerical  calculations  and  further  discretization  do  not  appear 
and  thus  their  effect  does  not  have  to  be  suppressed  with  further  filte¬ 
ring. 

The  object  to  be  reconstructed  is  assumed  to  be  within  a  circle  of  ra¬ 
dius  D.  The  device  can  be  considered  to  be  the  equivalent  of  a  raster 
of  parallel  rays,  which  are  rotated  about  the  axis.  Let  0  index  the 
angular  position  of  the  raster,  l  the  distance  from  the  axis  to  a 
parallel  ray,  and  let  the  location  of  a  point  inside  the  circle  be 
given  in  polar  coordinates  as  (r,4>).  (See  Fig.l).  Then  f  ( r ,  <t> )  is 
the  X-ray  absorption  coefficient  at  the  point  (r,l>).  Let  p(t,t)  be 
the  line  integral  over  f  taken  over  the  ray  indexed  by  (i.O). 

We  begin  with  the  Radon  inversion  formula  for  parallel  beams  as  quoted 
by  Herman  and  Naparstek  (1977),  equ.  (6). 


f ( r  ,4)  =  1 im 

e-o 


2  it  D 

/  / 

o  -D 


kc(r-0 


p(t,0  )df do 


(3.1) 


,  i,  i  -  -  **’***/*.>  •«***  *  -»w^w»6lfttaiDwn 
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where  kc(u)  ■  if  | u | >6  and  ke(u)aO  otherwise,  and 


4 1  *  r  cos (0-4>)  . 

We  first  consider  the  idealized  case  where  the  beams  are  infinitely 
narrow.  Then  the  data  { z ^ j >  consist  of  noise  contaminated  measurements 
of  p(4i ,Oj)  ,  that  is  , 


zij  =  P<Vsj)  +  6  i  j  i”-N , . . .  ,N 

j*l , . . .  ,M  . 

It  is  desired  to  estimate  f ( r , d> )  on  a  grid  of  points 
{rk’*j}’  from  this  data*  11  wil1  simplify  matters  if  we  let  <}> . * 2 tt j / M  , 
j*l ,2 , . . . ,M.  Without  loss  of  generality,  we  may  derive  our  formulas 
by  setting  $=0  since  the  formulas  for  <(>=ok  may  be  obtained  by 
relabeling  the  data. 

First  fix  0=0^ .  The  inner  integral  in  (3.1)  becomes 


4  ’  -  e 
1  im  / 
c-*o  -D 


p ( £  ,0j )  d£  + 


1 _ d_ 

F^T  d£ 


pU.O^dl 


(3.2) 


The  idea  is  as  follows.  One  first  obtains  a  cubic  smoothing  spline 
approximation,  call  it  p.(4,0j  to  p(4,0.).  p.(4,0.)  is  the 

.  .  .  ,  A  J  J  A  J 

m l  m in i  z e r  of 


J  J_N  (f  (  *,•  )~2i  j  >2  +  A/  (  f(4))2d£  (3.3) 

2 

in  Wgt-O.D).  p^(4,0j)  has  a  representation 

Px(4,Oj)  =  ak+bkl+ck42/2+dk43/3  ,  *e(4k,4k+1]  , 

where  a  k  * b  k  » c  k  and  dk  are  (for  fixed  X),  linear  functions  of  the 
data  z  •  j  ,  i=-N,...,N,  and  px  has  two  continuous  derivatives’in  4  . 

If  X  is  chosen  well,  then  under  circumstances  that  are  likely  to  be 
satisfied  here.^it  is  known  that  fyU.Oj)  is  a  good  estimate  of 
gj  p ( 4 , 0 j )  .  See  Craven  and  Wahba  (1979).  We  estimate  p(4,0j)  by 

dTPxu-0j)  x  W+V2  •  ^Wii  •  O-4) 

There  exist  coefficients  wk<  ,  wj^  and  wki  independent  of  the  data 
and  depending  on  X  and  (l^)  such  that 

1/  but  see  last  paragraph  below. 
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b‘  ■  &&  !«i 

H 

ck  ■  !tj 


d*  ■  !1i 


(3.5) 


These  coefficients  can  be  determined  and  stored  once  and  for  all 
(after  X  is  selected),  requiring  3(2N)(2N+1)  storage  locations. 

The  storage  requirements  can  be  reduced  at  the  cost  of  time  by  exploiting 
recursion  relations  between  the  {w^}  »  we  omit  the  details. 


By  substituting  (3.4)  into  (3.2)  the  inner  integral  can  be  evaluated 
analytically  and  the  limit  as  e-*-0  taken. 


First,  let  Z'=rcos0j  be  in  the  interior  of 
becomes 


Then  (3.2) 


N-l  4i  +  1  b.+c,£+dV 

l  !  ■J—— J“““ -  dZ 


i  =  -  N  Z  . 
\fm 


Z'-Z 


+  1  i  m  1  / 

€+0  V  Z 


Z'-e  b„ +cr,Z+dmZ 
m  m  m 


Z'-Z 


m+1  b+c_Z+d  r 
+  j  m  m  m 

Z'+e 


Z'-Z 


dz 


(3.6) 


*  Jr(0j)  .  say  . 

Upon  carrying  out  the  Indicated  integrations,  one  obtains 


w 


N-l 

l 

i  =  -N 


(b^+c.z'  +  dj z'  Z)Zn 


+  (bm+c„,t,+d  z1 
mm  m 


+  -X  {ci+2dir)Uurti)  +  \  XI  di(U'-'i)2-(f-,-ti+i)2)- 

(3.7) 

These  calculations  are  all  stable  except  possibly  for  the  two  terms  invol- 
Z ' - Z  Z'-Z 

ngZn(-; - ^_)  and  either  Zn(yn — r— 1 )  (if  Z'  is  near  em)  or 

1  '*'in  m 
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£  1  -  il 

ln{  (if  jt'  is  near  A  +1 ) .  We  give  the  details  for 

1  ‘znH-2 

A'  close  to  A  . 

m 

Let  %+rVh  and  1et 


A 1  =  A+6h 
m 


(3.8) 


where  0 < 5< 1 / 2  .  The  possibly  offending  terms  from  (3.7)  are 


O  A  '  *  A_  , 

£■  \  .  _  /  m- 1 


Since  the  cubic  smoothing  spline  has  continuous  first  and  second 
derivatives  at  Am,  we  always  have  the  relations 


(b  -b  ,)+(c  -c  ,)t  +(d  -d  ,)i2  =  0 
v  m  in- 1  ‘  '  m  m-1'  m  '  m  m-1'  m 


(3.10) 


(c  -c  ,)  +2(d  -d  , ) i  -  0 

v  m  m-1 '  '  m  m-1 1  m 

Substituting  (3.3)  and  (3.10)  into  (3.9)  gives  that  (3.9)  is  equal 
to 


2<Vdn-l)h2«2tn<T7*>  +  (bm+cmi'+V2)in(T^) 


(3.11) 


If  1/2<6<1,  a  similar  expression  may  be  obtained  by  summing  the 
mth  and  m+1  st  terms.  , 


Subtituting  (3.11)  into  (3.7)  one  obtains,  provided  0<6<l/2,  (and 
assuming  the  A^  are  equally  spaced) 


J 


r«V 


N-l 


(bs+( i+5)hc1+(i+6)2h2di )  An( 


6  +  (m- i )  ' 
5  +  (ni- i  -  1  ) 1 


N-l  ,  N-l  , 

+  h  l  (c,+2(i+6)d.l  +  l  l  4. [ 2 (m+6 ) h- ( 2 i - 1 )h‘] 
i=-N  1  1  ‘  i =  -N  1 


+2(dm'dm-l)h262ln(T+6) 


+  Vcm(m+5)h+Vm+6>2*n(Trf)  •  (3.12) 

Since  52An(-^~j)  and  An(j^|) -►  0  as  5-*-0  ,  this  expression  is 
computed  in  a  s  tra  i  gh  tf  orward  ma-nner  for  €  Q  <  6  <  1  /  2  ,  for  some  suitable 


eQ  ,  and  set  equal  to  0  if  0<5<eQ  •  A  similar  expression  is  obtained 
for  1 /2<6< 1  . 


Having  evaluated  J  (0.),  the  estimate  of  f(r,o)  is 
r  0 

i  M 

f(r,0)  =  tr  l  J_(0,)  .  .  (3.13) 

j  =  l  J 

Thus,  one  can  process  each  set  of  projection  data  (i.e.  the  data  for 
fixed  0j)  in  parallel.  For  each  j  one  collects  Zj  =  (z  n  j^  ’ 

computes  the  {b^},  (ck)  and  from  (3.5)  ,  J r ( 0 ^ )  from  (3.12)  or 

the  corresponding  expression  for  1/2<<5<1,  and  f(r,0)  from  (3.13). 

To  obtain  f(r,Op),  0pfO,  one  repeats  the  calculations  with  each  data 
set  Zj  relabeled  as  zj.p-  Note  that  the  coeffients  bk,ck  and  dk 
depend  only  on  z . .  They  can  be  computed  in  parallel  once  for  each  set 
of  projection  data  and  then  the  projection  data  discarded. 

The  regularization  parameter  here  i  s  x  (the  choice  of  eQ,  if  reasonable, 
is  secondary ) .  If  X  is  fixed  the  wfc.  of  (3.5)  can  be  stored. 

The  ultimate  choice  of  X  (or  several  values  of  X  to  provide  alterna¬ 
tive  pictures),  should,  of  course  be  chosen  by  examining  pictures  with 
competitive  X  for  their  medical  usefulness.  Since  X  controls  the 
smoothness-fidelity  tradeoff,  varying  X  is  likely  to  have  the  visual 
effect  of  bringing  the  picture  in  and  out  of  "focus".  A  too  large  X 
should  result  in  an  oversmoothed ,bl urred  picture  while  a  too  small  X 
should  result  in  an  overly  grainy  or  "streaky"  picture.  A  good  set  of 
candidate  X's  should  be  obtainable  at  the  design  stage  by  using  the 
method  of  ’generalized  cross  validation  (GCV),  on  data  from  typical  subjects 
with  the  parameters  (e.g.  number  of  photons,  number  (2N+1)  of  rays,  etc.) 
that  will  be  used  in  practice.  Transportabl e  code  is  available  for  doing 
this  (Merz  (1979),  Utreras  (1979),  Fleisher  {(1979)).  Given  X,  the  coef¬ 
ficients  wfc<  may  be  obtained  from  standard  spline  theory  (e.g.  Reinsch 
(1967)).  Numerical  results  on  the  estimation  of  the  derivative  from  noisy 
data  by  this  method  may  be  found  in  Craven  and  Wahba  (1979). 

We  now  consider  the  case  where  a  line  integral  approximation  to  the  data 
is  not  adequate.  Snopose  it  is  more  appropriate  tc  assume 


say.  Then  pxU,0j) 


zij  ■  /  • 

*i 

is  estimated  by  the  minimizer  of 
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,  N-l  *i  +1  ,  D  , 

4-  l  (  /  w^lJftDdl-z,  JZ  +  X  /  (f  "U))Zd?.  . 
c  i=-N  i.  J  -D 

If  w .(l)  is  taken  as  a  constant,  then  p^(t,0j)  has  a  representation 

PAU,0j)  =  ak+fekt+ckt2/2  +  cf  t3/3+ekt4/4  ,  He[*k,tk  +  1]  where  ak , 

\,?k,\  and  are  linear  functions  of  the  data  and  p^  has 

3  continuous  derivatives.  Expressions  for  the  ak>&k,ck,  cfk  and 

ek  can  be  obtained,  for  example  by  using  the  representation  for  splines 

given  in  Wahba  and  Wendelberger  (1979).  An  expression  for  J r ( G j ) 

is  obtained  by  adding  terms  corresponding  to  e.  to  (3.12).  w.(i) 

1 

can  also  be  modelled  as,  e.g.  a  trapezoid,  which  will  still  result 
in  a  piecewise  polynomial  representation  for  with  a  sufficient 

number  of  continuous  derivatives  .to  carry  out  a  similar  analysis. 

There  will  be  more  pieces  to  the  piecewise  polynomial,  however. 

4 .  On  the  relation  between  the  spline  transform  method  and  Tihonov 
requlari zation 

In  section  3.  we  have  discussed  a  new  method  for  the  numerical  inversion 
of  the  Radon  transform  which  essentially  consists  of  smoothing  the 
data  in  the  range  space  and  then  inverting  the  transform  analytically. 
Due  to  the  circular  symmetry  in  0,  if  one  obtains  p^(f.,0)  for 
0/Ol .... ,0^  ,  by,  e.g.  any  periodic  spline  interpolant  in  0  through 
P^(t,Ok)  >  k=l ,2 , . . . ,M ,  and  then  performs  the  integrations  of  (3.1) 
exactly,  the  result  will  be  the  same,  namely  (3.13). 

We  now  d.iscuss  the  relation  of  such  methods  to  Tihonov  regul  ar  i  za  t  i  on . 
Let  H  be  any  Hilbert  space,  let  Lt>teT  be  a  family  of  linearly 
independent  continuous  linear  functionals  on  H  and  define  the  opera¬ 
tor  K  by 

( Kf ) ( t )  *  g(t),  g(t)  =  Ltf  ,  teT  . 

Letting  X  be  the  range  of  K,  we  can  make  X  a  Hilbert  space  with 
the  norm 

INI*  =  inf  ||f|| 

*  feW 
Kf=g 

Now  consider  the  data  smoothing  problem:  Find  gcX  to  minimize 
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i  5  (9<t1)-zi)2  +X||g||J  .  (4.1) 

Letting  nt  be  the  representer  of  Lt,  it  can  be  shown  that  X  is 
a  reproducing  kernel  space  with  reproducing  kernel  Q(s ,t)=<n,»n*>u- 

S  t  n 

(See  Nashed  and  Wahba  (1974).  It  then  follows  that  the  minimizer 
gp  ^  of  (4.1)  is  gi ven  by 


9n,X(t)  =  (Q< t, .t) . Q(tn,t))(Q+nXI)"1 

where  Q  is  the  n*n  matrix  with  ijth  entry  Q( t • , t . ) =<nt 

J  t1 

Now  Q(t-,t)=Ltnt  =Qt  (t)  ,  say  ,  so  that 


z  ,  (4.2) 


> 


H  ’ 


Letting  L-=Lt 
minimizer  of  1 


satisfies 


Qt  =  Kn 

ti 

we  have  by  inspection  of  (2.3)  that  f 

1  I  (Lif-2i)2+X||f|g 
n  i  =  l 


the 


Now  nt  »  tcT  span  the  orthogonal  compliment  o(  the  null  space  of  K, 

since  <n  t » f  >  =  0,  tcT  Kf  =  0.  Thus  ffl  ^^(K)-1-,  so  that  fn  ^  is  the 

unique  element  of  minimal  norm  satisfying  Kf=g  and  so  (by  definition 

of  the  generalized  inverse  K  )  f  v=K  g  .  .  Thus  minimal  norm  smoo- 

n  n  i  a 

thing  in  the  range  space  (endowed  with  the  induced  norm),  with  an  exact 
inversion  is  equivalent  to  Tihonov  regularization  in  the  domain  space. 
Furthermore  E  |  !g  -  g  n  ^  1 1  ^  =E|(  K+  Kf  -  f  n  2  ,  where  the  expectation  is  taken 
over  the  and  convergence  obtains  as  n -►=>  under  general  conditions 

and  X  =  X ( n )  is  chosen  correctly.  See  Wahba  (  1  977  ). 


In  the  procedure  we  have  discussed,  smoothing  in  the  0  direction  is 
not  explicit  and  any  periodic  method  will  give  the  same  result.  Let 
PxU,0)  be  obtained  by,  say,  cubic  spline  interpolation  given  p^(Z,0k), 
k  =  l,2,...,M.  If  p  ^  { Jt ,  G )  were  the  minimizer  of,  say 


211  0  34d 


in™..)  '"'W-'u'  *  ‘  .1^^ 


d0d£ 


(4.3) 


in  an  appropriate  space  of  functions  periodic  in  0,  then  the  method 
being  proposed  would  be  exactly  equivalent  to  Tihonov  regularization. 


The  rainiraizer  of  (4.3)  and  px(i,0)  do  not  appear  to  be  exactly  the 
same  function,  however,  but  the  resulting  inversion  appears  to  be  close. 

As  far  as  the  choice  of  space  is  concerned,  this  method  is  appropriate 
for  p ( • , 0 ) e  w 2  , 

W*  =  { f : f , f '  continuous,  f “eL2  (-N ,N] }  . 

hk,.  'er,  it  is  more  natural  to  assume  p(  •  .OjeWg3  (f :  f  continuous, 

feLgl-N.N)}  as  follows:  Consider,  for  example,  head  sections  f(x,y) 

which  are  continuously  differentiable  functions  of  x  and  y  plus  a 

tumor  which  is  the  equivalent  of  adding  a  region  of,  say,  constant  higher 

density.  If  the  boundary  of  this  region  is  strictly  convex  and  "smooth" 

then  a  little  reflection  will  show  that  p(i,0)  is  a  continuous  function 

of  i  and  —  p  ( it ,  0 )  is  piecewise  continuous,  so  that  pf-.OjeW..  The 

preceeding  analysis  with  line  integrals  cannot  be  carried  out  to  obtain 

a  stable  computational  formula  because  the  derivative  of  the  linear  spline 

is  not  continuous.  However  a  similar  analysis  can  be  carried  out  with  a 

double  integral  over  a  0- i ncrement ,  or,  alternatively,  doing  spline 

smoothing  assuming  p (  ■  ,  * )  ewi { -N , N )  ®  (periodic).  Thiswill  appear 

^  ^  3 

se,.  -ately.  The  ability  of  —  px(i,0)  to  approximate  g-jrp(f.,0)  in 
the  w2  •->  has  yet  to  be  established  in  a  practical  sense  but  may  be 
quite  satisfactory  for  the  present  purpose  if  there  is  convergence. 


Appendix 


The  GCV  estimate  of  X  is  the  minimizer  of  V ( X )  given  by 

V(X)  -  I  jiUk(x)-zk)2/d  -  I  iIiak(X))2 

where  z k ( X ) f p > ^  ,  and  aR (X ) =  y|—  Lkfn  x  .  See  Wahba  (1979b)  and 
references  cited  there.  Letting  L  !^  =  L .  .  where  i  indexes  the  ray 
number,  and  j  indexes  the  rotational  position  of  the  detector 
(i  =  -N,...,N,  j  =  1 , . . . , H  In  the  notation  of  Section  3),  due  to  the 
rotational  symmetry  of  «the  device  one  should  have  a^ j(X)=a^(X)  indepen¬ 
dent  of  j ,  thus 


V  ( X ) 


-  I 

n  .  £■ 


i  =-N 


M  „  , 

[  (z.  ^Xj-z.JVd- 
j  =  l  J  J 


l  a  •  ( X ) ) 
i*-N 


(I  thank  F.  flatterer  for  this  observation). 


One  way  of  obtaining  the  denominator  is  to  compute  a.(x)  as  L... 

i  1  i  1  11 

6  ,  where  5  is  the  "picture"  when  the  input  is 

n ,  a  n ,  x 

z  =  (0, . . .  ,0,1  ,0, . . .  ,0)  with  1  in  the  ilth  position.  As  an  approxima¬ 
tion  to  V(\)  one  might  consider 

N  *  M  *  N 1 

V(x)=  /  2N  1  +1  )  m  1  £  ^zij^^’zij^  /  ^  *2N  ‘  +1  E  ai  ^  ^ 

k  =  _N.  l=]  Ve.  Vi  +i  k  =  .N.  ik 

where  { i k }  and  {j^}  are  either  representative  or  randomly  selected 
subsets  of  the  indices. 


We  note  that  the  use  of  GCV  is  appropriate  to  estimate  the  smoothing 
parameter  for  low  pass  filtering  methods  other  than  Tihonov  regulariza¬ 
tion,  see  Craven  and  Wahba  (1979),  Golub,  Heath  and  Wahba  (1979). 
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